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1 . Introduct ion 



This review is intended to summarize the developments in numerical weather 
prediction which have occurred during approximately the last three years* This 
roughly covers the period since the completion of two major texts by Haltiner 
(1971) and Phillips (1973). Earlier research will be referenced only when 
required for the development of a particular topic. An indication of the in- 
tense activity in this field can be seen in the papers presented at the Second 
American Meteorological Society Conference on Numerical Prediction in Monterey, 
California during 1-4 October 1973, for which the abstracts have been published 
in the July, 1973 Bulletin of the American Meteorological Society . 

Since the first application of numerical integration techniques to large 
scale weather prediction [Charney, Fjortoft, and von Neumann (1950)], these 
techniques have been applied to a large variety of atmospheric and oceanic pre- 
diction problems. These include studies of the general circulation of the 
atmosphere and the oceans, tropical storms, fronts, sea breezes, cumulus con- 
vection, atmospheric pollution, the planetary boundary layer and turbulence. 
This review will focus on the developments in large scale weather prediction, 
although many of the results have applications to other problems. The para- 
meterization of smaller scale phenomena will not be considered in this review 
because it deserves a separate comprehensive treatment. 

Numerical integration of the atmospheric equations as an initial value 
problem is the primary basis for the prediction of synoptic-scale disturbances 
for periods between 12 hours and perhaps five days and, in addition, to some 
extent for smaller scales and much longer periods. The sources of error in 
such predictions are a consequence of (a) gaps and errors in the data which 



4 



make up the initial state, (b) limitations in the objective analysis-initiali- 
zation schemes which are applied to the data, (c) truncation errors in numeri- 
cal integration schemes, (d) incomplete representation of the many complicated 
dynamical processes at work in the atmosphere and finally, limitations imposed 
by the predictability of the atmosphere* 

2. Objective Analysis and Initialization 

A standard method of objective analysis which has been in use for nearly 
15 years or so begins with an initial guess based on a prognosis of the mass 
field from the previous analysis or, perhaps simply a persistence forecast. 

The initial gridpoint values are then modified with actual observations 
weighted according to the distance between the observation and the gridpoint. 
Winds are obtained geostrophically or by means of a balance equation. With 
some care to insure vertical consistency and to avoid superadiabatic lapse 
rates, this type of objective analysis is quite adequate for initializing the 
vorticity (filtered) type of prediction models. However, the primitive equa- 
tion or P.E. models introduced operationally about a half dozen years ago, 
are much more sensitive to the initial conditions than are the filtered models 
from which gravity waves have been eliminated. As a consequence, the lack 
of proper balance between the pressure and wind fields can give rise to quite 
large, spurious inertial-gravity waves with periods ranging from a few time 
increments At to a day or more. Through the natural geostrophic adjustment 
process, which is inherent in the P.E. models, these spurious gravity modes 
are eventually dispersed and suppressed to acceptable magnitudes, perhaps 
comparable to those observed in nature. To avoid the large pressure fluctua- 
tions associated with these spurious gravity waves, greater care should be 
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taken to bring the objective analysis to a state where the wind and pressure 
fields are more nearly in proper balance before beginning the forecast, a 
procedure which is usually referred to as initialization * 

The first efforts in this direction simply followed the practice in fil- 
tered models by using winds from a stream function obtained by solution of 
the balance equation with the objectively-analyzed geopotential fields on 
constant pressure surfaces. Since the balance equation allows for curved 
motion to some extent, the resulting winds are in better balance with the 
pressure force than geostrophic winds. Nevertheless, as shown by Phillips 
(1960), an appropriate divergent wind component is necessary if the inertial- 
gravity waves are to be excluded or negligible. Consequently, the rotational 
winds obtained by solution of the conventional balance equation are not 
sufficient to prevent the unwanted gravity modes in a primitive equation 
model, although they are clearly better than geostrophic winds or unmodified 
observed winds. 

Recently, Sundqvist (1973), on the basis of limited tests, suggested that 
the solution of a balance equation with zero mass divergence (V • TT\V = 0) 
on sigma surfaces would lead to initial winds which give less gravitational 
noise than the winds from the conventional balance equation on pressure sur- 
faces followed by interpolation to sigma surfaces. 

In another approach, Temperaton (1973) verified that gravity oscillations 
were indeed present after initialization with the exact rotational part of 
the wind. The latter was obtained from an idealized experiment where ’Con- 
trol" data was first generated by first running a numerical prediction model 
for two days using the Euler-backward scheme to dampen the high frequency 
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modes leaving essentially only the slowly evolving meteorological modes. 
Then the geopotential field was treated as the observed field and a series 
of experiments were run for an additional day with different initializa- 
tion methods. In particular the rotational wind component was extracted 
from the control field by solving the Poisson equation, v 2 ^ = k • Vx V , 
for the stream function • Initializing with this wind resulted in 
gravity noise comparable to that resulting from winds obtained from solu- 
tion of the balance equation. 

The divergent component of the wind implied by the quasi-geos trophic 
theory can be obtained by solving the corresponding co- equation and then 
using the continuity equation to obtain the divergent wind potential X 
as follows: 

a- a* 

This is a rather lengthy task however, and is probably less effective and 
maybe as time consuming as simply using the model equations to achieve 
the desired balance. In fact, Houghton, Baumhefner and Washington (1971) 
made a study of initialization for the NCAR global model and found that 
inclusion of the vertical velocity from an co-equation and the corresponding 
divergent wind component in the initial winds did not significantly reduce 
the trauma of inertial gravity oscillations. 

At the National Meteorological Center (NMC) solution of the balance 
equation was abandoned several years ago in favor of extracting the rota- 
tional wind component directly from objectively analyzed winds. In addi- 
tion, the 12-hour forecast divergent wind component is utilized. The 
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resulting total initial wind field gives better agreement with the observa- 
tions without significantly increasing the noise level in the prediction 
model . 

Although they are usually treated as separate steps, objective analysis 
and initialization have the same goal, namely, to provide an accurate, 
properly balanced state from which the hydrodynamical equations can be 
numerically integrated forward in time. Obviously, the two steps are not 
dynamically independent and need not be separated, which is the view taken 
by a number of investigators. A specific example is Sasaki's (1958) ini- 
tialization by variational methods which incorporates dynamical constraints 
beyond the usual geostrophic wind law in the objective analysis. These 
constraints have included the use of a generalized wind equation, suppres- 
sion of high-frequency oscillations, as well as minimizing RMS differences 
between observations and a first guess field. 

Lewis (1972) applied Sasaki's technique to develop an operational analy- 
sis of wind and temperature for the global band 40°S to 60°N from the 
surface to 250 mb. Disregarding vertical advection, Lewis derived the 
following generalized thermal wind equation 

£ 55 ' - h * 55 - EVr B ? <2 > 

The variational formulation calls for the minimization of the integral 
I = fff [a(V - V) 2 + P(T - T) 2 + a(F x 2 + F y 2 )] (3) 

where a, P are precision moduli for the wind and temperature based on 
variances and a is a dynamical weight factor. Solution of the Euler 
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equations associated with the above integral provide the new wind and tem- 
perature fields. At the Monterey NWP meeting in October 1973, Sasaki 
applied a similar procedure to the problem of initialization by utilizing 
a dynamic constraint on the kinetic energy which reduced the gravity wave 
noise considerably. 

The variational principle has also been used by Flattery (1967) at NMC 
in making objective analyses by fitting Hough functions (the eigenfunctions 
of Laplace’s tidal equation) to the observed data by a least squares method. 
The results have been used in a global spectral model. 

Another common method of objective analysis interpolates between ob- 
servations to gridpoints by requiring the gridpoint values to satisfy a 
Poisson equation. Starting with an initial guess, the gridpoint values are 
then modified by the use of observed values which are treated as internal 
boundary points. This scheme was evaluated by Leary and Thompson (1973) 
for accuracy with known spherical harmonics. Wavenumber 2 was reasonably 
well represented with 87 $ of its amplitude squared appearing in the analysis 
field, while only 13$ of the input amplitude squared of wavenumber 12 sur- 
vived the analysis. Aside from demonstrating a rather severe deficiency 
of this analysis scheme, the study suggested that there is a less steep drop- 
off in the kinetic energy of the wind field, perhaps a "-2" rather than the 
"-3" power law obtained from spectra of such objectively analyzed data. 

Several other approaches have been utilized to suppress the inertial- 
gravity oscillations usually present in the early stages of a P.E. forecast. 
One procedure consists of integrating forward and backward about the initial 
time starting with the objectively analyzed data and periodically restoring 
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either the mass or the wind field at the initial time, usually the mass field 
in middle latitudes and perhaps the wind field in tropical latitudes. For 
this purpose it would appear advantageous to use an integration method which 
selectively damps high-frequency oscillations such as the Euler-backward 
scheme introduced by Matsuno (1966). As applied to a simple advective equation 



1-1 = 0 



(4) 



The Euler backward scheme will damp short waves about 10 $ with each time step, 
while long meteorological type waves will be negligibly damped. However, 
internal gravity waves of low frequency may also dampen rather slowly. Thus 
considerable computing time may be required by this two-step scheme to accom- 
plish the desired goal of suppressing the spurious gravitational noise, perhaps 
the equivalent of a 24- or 48-hour forecast which is operationally undesirable. 

The presence of the gravity modes is usually reflected in the horizontal 
divergence, the vertical velocity and the surface pressure tendency. Recog- 
nizing that the divergent part of the wind is intimately related to the 
propagation of the gravity oscillations, Talagrand (1972) suggested a special 
viscosity term to dampen the divergent wind component as follows: 
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Here D is the horizontal velocity divergence and Q , the relative vorticity. 
Now when the divergence and vorticity equations are formed, one obtains 
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It is clear that with V = 0 and (J. ^ 0 , divergence but not vorticity will 
be suppressed. McPherson (1973) of NMC evaluated this novel concept with an 
8-layer global primitive equation model using Shuman 1 s semi-momentum differ- 
encing. With (j = 0 , oscillations primarily in the form of external gravity 
modes appeared in middle latitudes as follows: 



AMPLITUDE 
5-10 meters 
20 - 25 meters 
20 - 60 meters 



PERIOD 
3-4 hours 
6-7 hours 
10 hours 



These waves were largely eliminated with a viscosity coefficient of |j = 2.5 

8 2 

x 10 m /sec in about 7 hours time. Several unfortunate side effects 

occurred ; however , there was abnormally high surface pressure at the end of 

12 hours in the vicinity of mountains together with a downstream trough. 

7 2 

Also values of \Jl in excess of 10 m /sec eliminated all precipitation. 
Clearly further evaluation is necessary. McPherson speculated that the 
slope of the sigma surfaces in the vicinity of mountains may be responsible. 
We would not find this surprising since two-dimensional diffusion on steep- 
sloping sigma surfaces near mountains may involve large vertical wind shear 
and temperature differencing, with undesirable consequences. Perhaps the 
diffusion technique would be more successful if carried out on quasi- 
horizontal surfaces, or in such a way as to avoid changing the selective 
vertical diffusion. 

At this point we would like to return to Temperaton's experiment on dy- 
namical initialization. He used the so-called "shallow water" equations in 
flux form with spherical coordinates. In previous experiments several years 
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ago Nitta and Hovermale (1967) had initialized by integrating forward N time 
steps, backward 2N steps, and then returning to the initial time with N steps 
forward; however, N was taken to be just 1. The Euler-backward scheme gave 
somewhat faster convergence than the leapfrog scheme. Mesinger (1972) recom- 
mended a larger amplitude NAt and partial restoration of mass more frequently. 
Winninghoff (1971) integrated a barotropic model 18 hours forward in time and 
back to the initial time with the 2-step Euler-backward scheme and obtained 
convergence to quite good initializing winds, but obviously at great expense 
in computer time (equivalent to a 72-hr forecast with the leapfrog scheme). 
Neither friction nor heating terms were present during initialization. Fric- 
tion and simulated heating after initialization gave somewhat poorer results 
than in the frictionless, adiabatic case. However, including these effects 
in the initialization procedure could lead to undesirable results depending on 
the mathematical formulation of these processes, which may not be reversible. 

Temperaton tried restoring mass only partially at the central time but 
was not successful. Using a time amplitude of 4At and 2/3 restoration of the 
mass field at every other time step gave quite rapid convergence with respect 
to the total RMS wind error. However, gravity wave activity, as measured by 
the RMS divergence, was not diminished. Analysis of the winds showed that 
the rotational wind component converged rapidly toward the correct value but 
the divergent wind component was more in error than with the case N = 1. He 
then tried another approach of integrating forward N time steps from the 
initial time, then N steps backwards from the initial time and finally averag- 
ing the end results as follows: 

u(t Q ) = j[u(t Q + NAt) + u(t Q - NAt)] (7) 
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The mass is now restored and the procedure repeated as many times as is needed 
to bring about convergence. The leapfrog method was used except for a first 
forward time step. The averaging tends to remove high-frequency modes but 
leaves low-frequency waves relatively unaffected. For example, the 2NAt fre- 
quency would be removed by such averaging while the 4NAt wave would be reduced 
by about one-half. Iteration cycles of 3 At, 6 At and 9At were tried with the 
best results obtained from the 6At (1 hour) case. The Euler-backward method 
was tried for starting in each direction but showed no improvement over a 
simple forward step. The principal result was that considerably faster conver- 
gence was achieved using the averaging technique than with the Euler backward 
scheme N = 1. The high frequency oscillations in the RMS divergence were 
markedly reduced to about one-tenth the amplitude occurring with the use of 
the exact rotational wind component without the averaging scheme. Nevertheless, 
there remained an easily recognizable oscillation with a period of about four 
hours. Elimination of this oscillation and further reduction of the noise 
was achieved by first adjusting the wind field while restoring the mass field 
at the initializing time and then repeating the forward and backward averaging 
procedure except that the mass field was averaged while the newly computed 
winds were restored after each cycle; thus both mass and winds were adjusted. 

In all, 20 cycles, equivalent to a 40-hour forecast, were completed, which is 
again too long for operational forecasting. 

Suppression of considerable gravity noise has been achieved at NMC by us- 
ing a running time average of the dependent variables in connection with the 
leapfrog scheme due to A. Robert (1966). Specifically, when a new value of 
a predicted variable, say A(t + At), has been calculated using the leapfrog 
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time differencing, a weighted average of the last three time values is cal- 
culated as follows: 

A* (t) = A (t ) + 0.5v[A*(t-At) - 2A(t) + A(t+At>] (8) 

The new value A*(t) is now used in the next leapfrog step to A(t + 2 At), 
that is 

A (t + 2At) = A* (t) + 2 At (§f) t+At (9) 

The weight factor v controls the amount of smoothing; V = 0.5 gives the 
familiar 1-2-1 averaging. The simpler averaging function with A(t - At) 
instead of A*(t-At) in Eq. (8) would completely remove waves of period 2 At 
and dampen a 4At wave by about one half, while much longer periods are rela- 
tively unaffected. Recall here that the spurious computational mode 
generated in the leapfrog scheme has very nearly a period of 2At if a corres- 
ponding physical mode is of a relatively lower frequency. As a consequence, 
NMC's averaging procedure is effective in removing the computational mode. 
Moreover, repetition of the averaging at every time step tends to damp some- 
what lower frequencies, including a considerable part of the spurious gravity 
noise generated through the early hours of a P.E. forecast. Asselin (1972) 
has computed the damping and phase shift characteristics as a function of fre- 
quency for a time filter of this form. Low frequency meteorological waves show 
negligible amplitude change and phase shift, but there is very effective 
damping of computational modes. 

In some recent initialization experiments with a five-level, global P.E. 
model, Haltiner and McCullough (1974) did not find any significant reduction 
of gravity noise with initial rotational winds from a balance equation solved 
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on C-surfaces as compared to similar winds from a p-surface balance equation. 
This may be somewhat surprising since from theoretical considerations, it 
might be expected that at least external gravity waves would be suppressed. 



Sundqvist assumed zero mass divergence, • TTV= 0 , where TT is surface 
pressure, to obtain the balance equation on sigma surfaces. It follows from 



if V • TTV = 0 , the surface pressure tendency will vanish initially which 
should tend to suppress external gravity waves. In any event, the Haltiner- 
McCullough initialization experiments with the a-balanced winds did not show 
significant improvement in damping spurious inertial-gravity waves. 

The Robert time filter (1966) is not only effective in eliminating the 
computational mode associated with the leapfrog scheme, but also gradually 
tends to dampen frequencies characteristic of the gravitational noise at the 
beginning of a P.E. forecast. Haltiner and McCullough combined the time 
filter with the averaging scheme of Temperaton to substantially reduce the 
noise in the equivalent computer time of a 12-hr forecast. The latter was 
accomplished by twice averaging winds from a 3-hr forecast and a 3-hr hind- 
cast and restoring the mass field to the initial values for each time. The 
winds thus obtained together with the original mass field constituted the 
initial conditions for prediction, which resulted in a substantial decrease 
in high frequency surface pressure oscillations compared to the immediate 
use of the balanced winds for initialization. 

The addition of the divergent wind component predicted by the model from 
a previous analysis to the present analysis time should further aid in damp- 
ing the spurious inertial-gravity noise generated by the introduction of new 
data. This step costs very little in computer time. 



the integrated continuity 




(TTv V) da , that 
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3. Data Assimilation 



The advent of the satellite infrared spectrometers (e.g., SIRS), which 
can sound the atmosphere from space, has provided a new and important 
source of data which is especially valuable over sparse-data areas. It 
did, however, further emphasize a problem of how to utilize most effectively 
upper air data which do not occur at or very near the regular synoptic obser- 
vation times, a problem already existing in connection with aircraft observa- 
tions. The procedure for incorporating off-time reports into an analysis- 
prediction scheme is referred to as four -dimens ional assimilation . Clearly 
the trauma of spurious inertial gravity waves associated with initialization 
must be eliminated if new observational data is to be more or less continuously 
incorporated into a prediction procedure. 

Some of the early idealized experiments dealing with the assimilation 
of SIRS data by Charney, et al (1969), Jastro and Halem (1970) and Williamson 
and Kasahara (1971) showed that the simple insertion of temperature data can 
eventually determine the wind field and vice versa. The process of inserting 
some variables into a numerical forecast while others are unchanged has been 
referred to as updating , which seems to be a rather restrictive definition in 
terms of operational forecasting. In any event, the degree to which one field, 
for example, winds, can be obtained by updating another, perhaps temperature, 
depends on the predictability of the model which, in turn, depends on natural 
instabilities, nonlinear exchanges, frictional dissipation, etc.. Using 
model data to stimulate the real world, then introducing errors and updating 
with, quote, "observed", data provides an indication of the best possible 
results that could be obtained by updating an operational prediction model 
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with real data. In actual operational forecasts there are other sources of 
errors including the representation of the physical processes and numerical 
truncation. Naturally, the better the model, the less the predictability 
error, but simple insertion of data is assuredly not the best way to update. 

It would be quite logical to modify other variables by static or dynamic 
balancing and refer to the whole procedure as updating or four-dimensional 
data assimilation. 

Mesinger (1972) suggested the following three names for methods of obtain- 
ing a proper balance between the mass and motion fields: 

(1) Static Balancing - use of a wind law such as the geostrophic 
relation or the balance equation. 

(2) Dynamic Balancing - going backward and forward about a central time, 
restoring mass or wind at t = 0 or letting both vary. Sasaki's variational 
method can also be considered as a combination of static and dynamic balancing. 

(3) Four -Dimensional Balancing in space and time - introduction of data 
three-dimens ionally periodically at discrete time intervals, including 
regular synoptic times. 

Hayden (1972) used a 2-layer P.E. model with Shuman's semi-momentum dif- 
ferencing scheme to run a series of experiments with the insertion of tempera- 
ture calculated from geopotentials which were in turn derived from the SIRS 
temperatures. Data were inserted six times during a 12-hr forecast. Tempera- 
ture tendencies over a period of one orbit were generally less than the 
observational error; consequently, as found by Talagrand (1971), too frequent 
updating may be deleterious. Hayden found that a static balancing, consist- 
ing of a geostrophic wind correction, Av, computed from the changes in 
geopotential A^ inferred from the temperature insertions aided the geostrophic 
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adjustment process and reduced the shock of updating. He also attempted 
dynamic balancing with N = 2 , i.e., integrating two steps forward, 4At 
backward and 2At forward to return to the initial time. 

Three criteria were used to measure the successful assimilation of data: 

(a) the updating must not unduly shock the model which he inferred from a 
measure of the divergent wind component; (b) does the model remember the data 
inserted (this was tested by reversing the forecast after 12 hours and noting 
whether the hindcast temperatures were nearer to the inserted values than 
were the original temperatures that existed prior to the updating.); (c) lastly 
does the updating result in better mass and wind distributions when compared 
to the NWS analyzed fields? For this purpose the 12-hour forecast was 
followed by a 12-hour hindcast. Then the difference between the cycled data 
and the initial data were correlated with the difference between the NWS 
analysis and the initial state; a positive correlation indicates success. 

The following conclusions were reached: 

a. Even with poor initial conditions, temperature can be assimilated 
without shock if some balancing is performed to aid the geostrophic adjustment 
process. Dynamic balancing is sufficient but static balancing in the form of 
a simple geostrophic wind correction speeds up the assimilation process con- 
siderably, at least outside the tropics. 

b. Hayden anticipates that under operational conditions, where the model 
state is maintained similar to the observed state, four-dimensional assimila- 
tion can be accomplished without time-consuming dynamic initialization. 

c. Four-dimensional data assimilation is evidently more effective than 
regular objective analysis of off time reports that have been updated to 
regular synoptic times by Lagrangian advection. The exact details of the 
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surface pressure appear to be relatively unimportant to the effectiveness 
of the four- dimensional assimilation. However, the SIRS-B data by them- 
selves are not capable of defining the circulation because the data density 
is not sufficient for objective analysis, nor with the model, of producing 
even an approximate surface pressure field. 

Bengtsson and Gustavsson (1971) had previously found also that analysis 
prior to insertion of data, say from a satellite, leads to a more rapid 
reduction of error. Talagrand and Miyakoda (1971) showed that a synthesis 
technique of averaging forecasts made from objective analyses made at 
different times can reduce the random errors of measurement and analysis, 
sort of a Monte Carlo approach. They did some studies of inserting data 
into a running forecast when and where the data were available. If the 
difference between the predicted values and the observations is less than 
or equal to the observational error, don* t ins ert it; there* s no point to 
shocking the system and uselessly creating spurious inertial-gravity waves. 

In summary, experimentation thus far suggests that one or several of 
the following steps are useful during initialization and four -dimensional 
data assimilation. 

Firstly, at the regular 0000Z and 1200Z synoptic times: 
a. Perform a preliminary objective analysis and obtain a stream func- 
tion by solution of the balance equation, perhaps on the surfaces of the 
vertical coordinate, or as an alternative obtain the rotational wind directly 
from objectively analyzed winds. A combination of both may be desirable, 
the former in middle latitudes and the latter in the tropics. 

In the tropics where the observational errors in the pressure field may 
be comparable to the pressure variations associated with synoptic disturbances 
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(except for tropical storms), it is generally unwise to attempt to obtain 
the stream function from the pressure field by solving the balance equation, 
which, in fact, is singular at the equator. A recommended procedure is to 
calculate the vorticity directly from the observed wind field and then obtain 
the stream function by solving the Poisson equation = Q . Finally, the 

geopotential field is calculated from the stream function with some form of 
the balance equation. Saha and Suryanarayana (1971) made a series of calcu- 
lations of the geopotential in this manner from the quasi-geostrophic rela- 
tion, the linear balance equation, the balance equation and the so-called 
vorticity form of the balance equation which are, respectively, 

V^(j) = fV^ \jj 
v 2 ^ = V • (fVi|f) 

( 10 ) 

+ v • w 

v 2 ^ = + k • vqxv - vV 2 /2) 

Observed winds were used to evaluate the vorticity, Q = k • VxV , and \|f was 
obtained from V 2 !)/ = J . The geopotential fields obtained from the last 
three forms were very similar and compared favorably with the analyzed geo- 
potential fields at the 850, 700, 500 and 300-mb levels. However, the last 
equation, the vorticity form, gave the least RMS error. 

b. To supplement the rotational wind component, the divergent wind com- 

ponent predicted from a previous analysis may be used. The cost in computer 
time amounts to the solution of a Poisson equation, V^X = V • V , followed 
by the calculation V = . 

c. Next apply some dynamic balancing such as Temperaton f s technique. 
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Step a and perhaps b may be replaced by a variational method of analysis in- 
cluding some dynamical constraints. 

Between synoptic times , insert new data where available at intervals of 
not less than three hours with an objective analysis in a limited region of 
the reports including first static balancing to adjust the wind or the mass 
fields partially to one other. This may be followed by dynamic, balancing to 
reduce further the gravitational noise if the differencing scheme and fric- 
tional terms are insufficient for this purpose. 

4. Integration Methods 

a. Semi- implicit schemes 

Although the momentum or primitive equations are simpler and involve 
fewer approximations than the filtered equations, the presence of gravity 
waves requires a much smaller time step to avoid computational instability 
with explicit integration schemes. Otherwise the small step is of little 
advantage or even harmful insofar as meteorological waves are concerned, and 
it is certainly expensive in terms of computer time. As a consequence, there 
has been considerable effort to circumvent the stability requirement. In the 
Soviet Union, Marchuk (1965) introduced a differencing scheme which treated 
the gravity modes implicitly and the low frequency meteorological modes 
explicitly, thus permitting a much larger time step. Kwizak and Robert (1971) 
successfully applied a semi- implicit differencing method similar to one 
suggested by Kurihara (1965) to a barotropic 500-mb forecast. To illustrate 
the technique we adopt the notation of Elvins and Sundstrom (1973) applied 
to the shallow water equations, which in differential form are 
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( 11 ) 



u = - (b - u u - v u + fv 
t T x x y 

v = - 4) - u v - v v - fu 
t T y x y 

(j) = - (|) (u + v ) - (u{) + vd) ) 

Y t r x y r x T y 

The pressure gradient terms in the momentum equations and the divergence term 
in the continuity equation, which together primarily govern the propagation 
of gravity waves, are evaluated semi- implicitly while the remaining terms 
are evaluated explicitly. 

In the difference equations u n+ ^ and v n+ ^ in the continuity equation 

are replaced by substitution from two momentum equations obtaining an elliptic 

. . in~hl 

equation m (p 

[l-At 2 f (D qx 2 + D qx 2 )] (f +1 - (j 1 ) = F (x,y,t n ,t n-1 ) (12) 

Here F is composed of terms at times t^ and t ^ • After solving this equa- 
tion for by one of the usual methods, e.g., relaxation, u 11 ^ and 

n+ i 

v may be calculated directly from the momentum equations. As a consequence 

N 

of the implicit character of the difference equations with respect to gravity 
waves, a much larger time step is permitted without encountering linear com- 
putational instability. Of course, extra time is taken at each step to 
solve the elliptic equation, but overall, computer time is reduced by a factor 
of 3 to 4. 

Elvius and Sundstrom (1973) suggested an efficient differencing 
system which is staggered in both space and time. This scheme not only 
permits a larger time step but also reduces the phase speed error of the 
low frequency or "meteorological" mode. They also developed suitable 
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boundary conditions for use with a fine mesh model, which is undergoing fur- 



ther tests with realistic initial conditions. Treating the coriolis terms 
implicitly permits a slightly larger time step. 

Baroclinic models are more complicated, but the semi-implicit technique 
is applicable in a similar manner. Gerrity, McPherson and Scolnik (1973) 
have developed the semi- implicit difference equations using Shuman's semi- 
momentum differencing technique for the NMC 6-layer primitive equation model. 
The model has run stably for four days; however, it is not yet operational. 

When the implicit scheme is used, the phase velocities of gravitational 
oscillations are greatly reduced, hence the geostrophic adjustment process 
can be retarded, as will any initialization procedure for damping spurious 
inertial-gravity waves. Experiments by McPherson and Kistler (1973) verified 
the delayed damping of the gravity waves; however, because of the larger 
time step, dynamic initialization can be accomplished with less computer 
time. The net result was still a gain of 2 to 1 over the explicit integra- 
tion during an initialization procedure. 

5. Direct Methods for Helmholtz and Poisson Equations 

The need to solve Helmholtz-type equations in connection with the semi- 
implicit methods has stimulated interest in the more recent "direct" methods 
of solving Poisson and Helmholtz equations. Leslie and McAvaney (1973) have 
compared the speed and accuracy of a number of methods for solving equations 
of the form: 



where the Helmholtz coefficient a and the forcing function f are known. 




(13) 
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When approximated by the usual finite difference analogues, the foregoing 
equation is expressible as a system of linear equations with a block tri- 
diagonal matrix of coefficients. In the past, iterative methods such as 
the Liebmann successive over-relaxation method ( SOR ) have been widely used 
because of their simplicity. However, faster direct methods developed in 
recent years are replacing the relaxation methods, at least when a rectangular 
region is involved with simple boundary conditions. The direct methods may 
be divided into four categories as follows: 

a. Block Method which uses the fact that A is block tridiagonal. 

b. Cyclic Reduction Method (DCR) which reduces the dimensions of the 

matrix to be solved in a recursive manner. This method is restricted 
to certain numbers of interior points, such that with an MXN grid 
either M or N must equal 2 - 1 . 

c. Matrix Reduction which through coordinate transformation reduces 
the problem to a simpler tridiagonal form that has an easily com- 
puted solution. The dimension reduction method ( DRM ) is relatively 
simple when a fast Fourier transform is available. 

d. Finally, in certain circumstances the fast Fourier transform can be 
applied directly. 

Table 1 shows some comparative solution times and accuracies (RE) of 
the DCR, DRM and SOR methods applied to Poisson and Helmholtz equations on 
a 65 x 65 grid. 
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TABLE 1 



(from Leslie and McAvaney) 
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It is clear that with respect to Poisson type equations, the direct methods 
are far superior. To solve the Helmholtz equations by direct methods, the 
Helmholtz coefficient has to be a constant. On the other hand, the itera- 
tive (SOR) methods permit a variable Helmholtz coefficient and give rapid 
convergence when that coefficient is large. This technique can then compare 

favorably with the direct methods, particularly if a little less accuracy is 

-4 

acceptable, say, RE ~ 10 . Also the SOR methods, which are most readily 

generalized to irregular domains and mixed boundary conditions, are simple 
to program and require minimal storage in comparison to direct methods. 

But it is hoped that the lack of flexibility in the direct methods will be 
overcome in view of the importance of solving Helmholtz equations in connec- 
tion with semi- implicit methods which can be real time savers with respect 
to computing. 

6. Global Grids 

As computing capability has improved in the last several years efforts 
have increased toward operational global forecasting and also in fine mesh 
models. The absence of lateral boundaries in a global model is an important 
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advantage since such unnatural barriers cause errors in hemispheric models 
that eventually propagate from the fictitious boundaries in the tropics into 
middle latitudes. Numerical forecasting in the tropics, which admittedly 
has severe limitations at present, would be quite hopeless if artificial 
boundary conditions are imposed in low latitudes. On the other hand, a 
tropical band with middle latitude walls is not satisfactory either because 
these boundaries are far too active. So the only alternative appears to be 
global models despite the vast areas of little or no conventional data in 
the Southern Hemisphere, albeit the weather satellites are helping to over- 
come the data problem. 

A significant difficulty with global models is the lack of a suitable 
plane projection which does not seriously distort some areas. The most 
natural approach is a latitude-longitude grid; however, the convergence of 
the meridians poses a knotty problem, for as the distance between equal longi- 
tude spacing shrinks toward the poles, a shorter time step is needed to main- 
tain linear computational stability. A time step short enough to maintain 
stability in polar regions is exceedingly wasteful in low latitudes. Various 
techniques have been used to overcome this difficulty with varying degrees 
of success. The most common approach of late has been to filter out or 
dampen the short waves that would lead to instability near the poles so 
that a relatively large time step can be used throughout, rather than simply 
decreasing the time step with increasing latitude. In the Arakawa (1972) - 
Mintz model, the procedure consists of temporarily Fourier analyzing fields 
which will be differenced with respect to longitude and then modifying such 
Fourier amplitude so that the CFL stability criterion can be satisfied 
without shortening the time step. This principle can be illustrated with a 
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simple advection equation for which the stability criterion for wave number 
k and maximum wave speed c is typically of the form 



sin kd < e (14) 

j 

where d. = a cos cp. M is the grid distance at latitude cp. and e is a 
J j J 

constant perhaps unity. By reducing the amplitude of each wave component 

of the longitudinal gradient at each latitude by a factor S., the criterion 

J k 

becomes 



(S \ 

v jk aAA. cos cp. 



sin kd ^ C 



(15) 



It is clear that for a fixed and At , computational stability can be 

maintained by decreasing S as the latitude and wave number increase. 
This type of procedure is applied to all longitudinal derivatives in the 
terms involving gravity wave propagation. 

Vanderman (1970) used running averages of the tendencies to filter high 
frequencies components; however, this can cause computational instability. 
At the NOAA Geophysical Fluid Dynamics Laboratory Holloway, Spelman and 
Manabe (1973) applied space filtering to all time integrated variables at 
each time step. The filtering limits the east-west wavelength at all lati- 
tudes to the distance of two gridlengths at the equator. The minimum wave- 
length is given by ^ m £ n = 4-TTa N ^ , where N is the number of gridpoints 
around a latitude circle and a is the radius of the earth. The maximum 
number of waves at latitude 0 is 



2TTa cos 8 

L . 
mm 



(N/2) cos 0 



(16) 
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which is rounded to the nearest integer • This is accomplished by a Fourier 
analysis followed by synthesis with only the desired waves at each latitude 
circle included. The procedure has no significant effect on quadratic con- 
servation properties of the differencing schemes. The components of vector 
variables are first transformed into polar stereographic coordinates to 
avoid the problems associated with averaging vectors from widely separated 
longitudes where the unit vectors differ substantially in direction as dis- 
cussed by Shuman (1970). 

When tested on barotropic and baroclinic models the foregoing procedure 
proved to be superior to the Kurihara grid which has a poleward decrease in 
the number of gridpoints per latitude circle in such a way as to keep the 
distance between gridpoints from decreasing appreciably. The problem of 
spurious anticyclogenesis at the poles associated with a Kurihara grid was 
diminished. It was concluded that fast Fourier methods and new parallel com- 
puters can provide the necessary speed to handle the extra gridpoints in 
polar regions when the simpler differencing schemes are used with spherical 
coordinates . 

Some tests were also conducted by Williamson and Browning (1973) with 
global grids which verified this conclusion. They found with a grid that 
is uniform in a curvilinear coordinate system, the accuracy of approxima- 
tions involving curvilinear velocities is less near the singularities. In 
order to avoid the small time step associated with a uniform grid, they 
tried the method of skipping points near the poles to maintain a more nearly 
uniform distance between gridpoints, but the skipped grid resulted in large 
errors. More accurate interpolations did not help this matter. However, 
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by applying the Fourier technique to remove short wavelengths the errors 
were comparable to a uniform grid requiring a much shorter time step. 



7. Fourth-Order Differencing 

In another aspect of computat ion-time vs accuracy, it has been shown 
that a greater reduction in phase error can be achieved per unit of extra 
computer time by selective use of fourth-order space differencing than by 
reducing the grid size [see, for example, Williams (1972)]. This may be 
illustrated with the simple equation (4). Suppose At is chosen to main- 
tain linear computational stability in a P.E. model which permits a fast 

external gravity wave with a phase speed C of perhaps 300 m/sec or more 

S 

and C At /Ax ^ 1 • For the slower meteorological waves, with phase speed, 
8 

say C , the ratio C At/Ax will be much less than one, perhaps a tenth, 
m m 

or so. Some computations were made for a centered difference form of (4) 
with both second- and fourth-order space differences and second-order time 
differencing (leapfrog). As an illustration of the improvement in phase 
speed accuracy, assume At/Ax = 0.2 . Then the ratios of the finite 
difference phase speeds to the true speed for several wavelengths are as 
follows : 
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2nd order 


0.64 


0.84 


0.91 


0.94 
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4th order 


0.86 


0.97 


0.99 


1.00 


1.00 



Although this is a much simplified illustration, improvements in phase 
speed accuracy of 5 to 20^ or more can be expected with fourth-order space 
differences. Moreover, the latter need only be applied to the terms govern- 
ing the propagation of meteorological waves and not the terms involving 
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gravity wave propagation* It should be mentioned that the fourth-order 
approximation required a somewhat smaller time step, perhaps 20 to 25^, 
to maintain linear computational stability [see Haltiner and Williams (1973)]. 
On the other hand, halving the grid distance in the simple one- dimensional 
model described above would quadruple the computation time, yet the re- 
sulting improvement in phase speeds would be roughly the same as going 
from second to fourth order, as may be seen by comparing the 4Ax and 8Ax 
ratios above (0.64 and 0.91), or the 6Ax and 12Ax values, 0.84 and 0.96. 

8. Mesh Models 

Because of the enormous range of scales of atmospheric phenomena and 
limitations of even the most modern computers, it is obviously impossible 
to model numerically all phenomena on a single mesh of uniformly spaced 
gridpoints. Scales which are too small to be represented with a specified 
length must be parameterized in terms of the large-scale variables if 
their influence is to be felt. On the other hand, influences from outside 
the region of integration can be imposed through the boundary conditions, 
but both procedures may be inadequate at times. It is noteworthy that 
some important mesoscale phenomena may be infrequent in occurrence and be 
predictable for only a short range of time. Consequently, it is desirable 
to superimpose a fine mesh grid (FMG) on a coarse mesh (CMG) covering a 
much larger area, perhaps a hemisphere or the entire globe. Quite a few 
meteorological organizations, both foreign and domestic, are currently 
carrying out numerical integrations on fine-mesh, limited area grids. 

One of the most critical problems in dealing with limited area fore- 
casts, including the superposition of different grids, is the treatment 
of the boundary conditions. This problem is not really new and had to be 
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faced in the first numerical prediction experiments by Charney, Fjortoft 
and von Neumann (1950). They concluded correctly that the optimal pro- 
cedure was to specify precisely as many boundary values as required by 
the corresponding linearized equation to have a well-behaved solution. 
Additional values needed for the finite difference equations should be 
computed by extrapolations from interior values. Apparently their method 
of extrapolation proved to be unstable, as later shown by Platzman (1954), 
and they simply specified all values on the boundary. This gave stable 
results which, however, were less accurate and more stringent than neces- 
sary. Having to maintain constant values along inflow boundaries can lead 
to large errors propagating into the forecast region; however, this is 
certainly less of a problem when fine-mesh boundaries are permitted to 
change periodically through the use of coarse-mesh forecasts. Neverthe- 
less, the latter situation is in a sense a special case of the limited 
area forecast, for although the fine-mesh boundary values are no longer 
constant, it is necessary to obtain them by interpolation in space and 
time from coarse-grid values. The objective then is to do this in such 
a way as to avoid any instability at the boundaries and to obtain the most 
accurate forecast possible. 

The proper procedure for a barotropic primitive equation model, 
proposed by Elvius and Sundstrora (1973), following Charney (1962), is to 
prescribe at all boundary points the quantity > which corres- 

ponds to the outgoing characteristic. This permits gravity waves to leave 
the region rather than be reflected. In addition, the tangential compo- 
nent of velocity is prescribed at inflow points. Over-specification of 
boundary values is less accurate and may lead to parasitic waves and 
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perhaps instabilities. For a baroclinic system the situation is much more 
complicated. Sundstrom (1973) recommends an approximation which specifies 
the tangential velocity at inflow points as above while at outflow points 
the tangential velocity is extrapolated from the interior. Next a com- 
bination of normal velocity component and potential temperature corres- 
ponding to the inward external gravity wave is specified on all boundaries, 
and the combination corresponding to the first outward gravity wave is 
computed by extrapolation from interior values in such a way as to avoid 
instability. This would not allow the internal gravity waves to leave the 
region. 

Returning now to actual experiments with nesting of grids, the early 
efforts of superimposing a fine mesh on a coarse mesh consisted of inte- 
grating the coarse mesh model, for say 24 hours, and saving the "history 
tape" of data at every time step. Then the boundary values for the fine 
mesh model are interpolated in space and time from the CMG predictions, 
which is an overspecification and less accurate but nevertheless is compu- 
tationally stable [Elvius and Sundstrom (1973)]. Here the FMG predictions 
are influenced by the large scales depicted on the CMG, but not vice versa. 
Nevertheless, the fine mesh gives better resolution of physical features 
and permits smaller scale disturbances to develop. Some examples of this 
approach are the experiments of Hill (1968), Wang and Halpem (1970). 

A successful meteorological experiment which permitted interaction 
between the CMG and FMG through simultaneous integrations was carried out 
by Harrison and Elsberry (1972). By utilizing a movable fine mesh centered 
over a traveling one- dimensional gravity wave the FMG produced forecasts 
comparable to those obtained by a fine mesh everywhere. The boundaries 
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of the adjacent CMG and FMG overlapped as illustrated below. 
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In their scheme the FMG integration determines the values at points 6 and 8; 
then these values are used for subsequent predictions for points 4 and 10 in 
the CMG. Similarly, the FMG predictions for points 5 and 9 utilize the 
values at points 4 and 10 determined by the CMG integrations. Thus, for 
example, with a simple advective equation chi/dt + U chi/cbc = 0 , the inter- 
face equations are 



5T = " 2d (u 6 " u 3 } ’ and (17) 

Bu 5 U 

5tT = " T (u 6 ~ V (18) 

The computed changes for the larger time steps at the CMG points which form 
the boundaries of the FMG (i.e. points 4 and 10) are proportioned equally 
to supply the intermediate temporal values of u^ which are needed at the 
smaller time intervals to integrate the adjacent FMG points (i.e. points 5 
and 9). Similarly, linear spatial interpolation is performed to obtain 
values between CMG points as needed for the FMG integration. As a conse- 
quence of this procedure the FMG integration influenced the CMG, and very 
importantly the integration procedure is stable and created no significant 
noise or near-discontinuities at boundaries. 
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Applying this procedure to an idealized two-dimensional tropical disturb- 
ance and keeping the fine mesh grid centered over the disturbance gave results 
equivalent to a fine mesh everywhere in terms of forecast intensity and also 
for energy fluxes across the boundaries dividing the grids. However, some 
’'noise" did develop at the boundary. This was suppressed by smoothing across 
the boundary and the addition of a small diffusion term. 

Phillips and Shukla (1973) considered the two strategies of one-way inter - 
action using the history tape and the two-way interaction such as the one just 
described. By a heuristic argument based on the method of characteristics, 
they inferred that the two-way interaction procedure would give a more faith- 
ful reproduction of the proper transmission of information into and out of 
the fine-mesh region. Some numerical tests with the shallow-water equations 
showed that the two-way strategy did indeed lead to less error. They also 
found that the Lax Wendroff two-step scheme gave a somewhat larger error at 
12 hours than did the leapfrog scheme, but the reverse was true at 24 and 48 
hours . 

Ookochi (1972) combined a fine mesh with a coarse mesh for the integration 
of barotropic primitive equations in flux form on a staggered grid. The re- 
sults were essentially a composite of complete fine-mesh and coarse-mesh 
integrations with no significant noise at the boundaries. The principal 
integral properties involving mass, total energy, etc. were well conserved 
during the 96-hour experiment. 

Harrison (1973) describes some further experiments with systems of two 
and three nested grids for the simulation of a tropical storm by integration 
of the primitive equation on a four-level model. His calculations demonstrated 
the feasibility and advantages of nested grids in savings of computer time. 
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Presumably better phase speeds would be achieved as a consequence of the 
fine mesh because of smaller truncation errors. Note, however, that any 
wave, particularly those that are poorly represented on the coarse grid, 
will change phase speed when passing through the interface into the fine 
mesh, and again later when it leaves the fine mesh. As a consequence, 
erroneous interaction can occur with that part that remains in the coarse 
net . 



9. Vertical Coordinates 

The most commonly used vertical coordinate in primitive equation pre- 
diction models is the sigma coordinate (a) which leads to a = 1 and 
da/dt =0 on the lower boundary. This simplification of the lower boundary 
conditions is accompanied by a more complicated expression for the pressure 
gradient force. 

The use of potential temperature as a vertical coordinate [Eliassen 
(1962)] has received new attention in the last few years. If there is no 
heating the isentropic surfaces will move as material surfaces. This feature 
is very helpful for resolving frontal zones and sharp upper level jets. 

A possible disadvantage of the isentropic coordinates is the fact that the 
potential temperature shows a large variation along the lower boundary. 
Primitive equation integrations have been carried out by Eliassen and 
Raustein (1967), Gall (1972) and Shapiro (1973) with isentropic coordinates. 
They encountered no major difficulties, and they obtained good simulations 
of the jet stream structure. Gall (1972) included latent heat in his study. 
Bleck (1973 a,b) made forecasts with real data with the quasi-geos trophic 
form of the isentropic potential vorticity equation. The forecasts showed 
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skill in predicting cyclone development. The isentropic models show con- 
siderable promise for limited area models which treat the smaller scale 
synoptic features. These models have not been tested in global formulations 
or for longer period forecasts. 

10. Forecasting Skill 

With the greater sophistication in numerical modelling and the large 
advances in computer technology, it may be rightly asked whether there have 
been concomitant improvements in forecasting skill. To answer this question 
we may turn to verification statistics published by the National Weather 
Service which are probably indicative of other groups as well. Figure 1 
from Shuman (1972) shows the S^ scores, which are approximately a measure 
of the normalized RMS vector error of pressure gradient, for 30-hour, sea- 
level (upper scale) and 500-mb (lower scale) forecasts for North America 
from 1948 to 1971. (Shuman states that in terms of practical skill, scores 
of 0.30 at sea level and 0.20 at 500 mb are nearly perfect while scores of 
0.80 at sea level and 0.70 at 500 mb are considered essentially worthless.) 
There is clearly a general downward trend over the years indicating increas- 
ing skill. The improvement in the latter half of the 50' s and early 60* s 
is ascribed to the introduction and continuing improvements in the barotropic 
numerical 500-mb forecasts, while at sea level there was increasing skill 
on the part of prognosticators in the use of the 500-mb forecasts. 

The first successful baroclinic model, a three-level filtered version, 
became operational in 1962, which is not reflected clearly in the graphs 
because the years 64-66 proved to be a more difficult period to forecast and 
barotropic forecasts suffered as well. Nevertheless, a new plateau of skill 
was established and there was a general downward trend of scores. The 
first surface baroclinic model was a simplified graphical type due to 
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Dr. Richard Reed which went into operation in 1962 and drew from the in- 
dependently-made numerical 500-mb predictions. The continuing improvement 
in surface predictions from 1962 to 1966 was largely a result of Reed's 
model . 

The NMC 6- layer primitive equation has been under continuous develop- 
ment since 1966; it is the first numerical model at NMC to produce a better 
sea-level prediction than manual predictions without NWP guidance. However, 
with the numerical product in hand, man can improve the prognosis by about 
five points on the average in terms of scores. The net result has been 
a continued improvement in the skill scores. 

Figure 2 from Cooley and Derovin (1972) shows the NMC average 36-hour, 
500-mb scores. The light shading shows the human improvement. Figure 3 
shows the 30-hour surface scores. 

Considering the marked increase in accuracy of surface and 500-mb prog- 
noses the corresponding improvement in forecasts of weather elements is 
somewhat disappointing. Precipitation-no precipitation forecasts for the 
periods 12-24 and 24-36 hours improved only about 5$ between 1959 and 1970. 
There was somewhat greater improvement in temperature forecasts. For ex- 
ample, the annual number of maximum temperature forecast errors equal to 
10°F or greater at Salt Lake City decreased from 60 in 1949 to 22 in 1971. 

Sanders (1973) recently reported on six years of forecasting tempera- 
ture and precipitation by staff and students of Massachusetts Institute 
of Technology for the local observation site at Logan Airport. Despite 
continuous improvement in predicted synoptic patterns at the surface and 
aloft, there was no increase of skill in temperature and precipitation 
forecasts as measured by the incremental accuracy over the climatological 
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mean. As may be expected, skill decreased with increasing length of fore- 
cast and more rapidly with precipitation than with temperature. Minimum 
skill occurred in summer, probably due to the greater influence of meso- 
scale phenomena. Sanders also indicated that the lack of improvement in 
forecasting temperature and precipitation and even a slight downward 
trend in the latter is also shared by NMC forecasters. He suggests that 
the reason for the first day forecasting difficulties may lie in the fact 
that the benefits of short-range synoptic-scale forecasts of the mass 
field have been maximized and that the source of the errors in precipita- 
tion and temperature forecasts are primarily a result of mesoscale 
phenomena such as fronts, land-sea circulations, convection, urban in- 
fluences, etc. According to Sanders* results skill in precipitation fore- 
casts drops to 10^ in 2.5 days and in four days for temperature. He 
suggests that something of a breakthrough in synoptic forecasting is 
needed to improve significantly prediction beyond the first day. 

Despite some variations in the trends of forecasting skills, it is safe 
to conclude that, overall, forecasting ability has shown definite improve- 
ment since the advent of numerical weather forecasts. Moreover, further 
improvements can be expected in the latter half of the 70* s as the current 
research efforts in numerical techniques, simulation of the physical pro- 
cesses, and initialization techniques become operational, along with better 
satellite data and smaller mesh sizes. 

We may now ask the question is there no limit to the ultimate range 
and accuracy of weather forecasts if one is willing to spend enough money 
to provide the necessary data and enough computing power for the calcula- 
tions? To answer this question the predictability of the atmosphere must 
be considered. 
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11. Atmospheric Predictability 

Lorenz (1965) pointed out that perfect prediction can never be expected 
because (a) the governing laws are not perfectly known but only approxima- 
tions, (b) nor is the system of equations strictly deterministic and finally 
(c) the initial state is not perfectly known. Moreover, he shows that even 
for a determinate system of equations, such as one normally used in numeri- 
cal weather prediction, separate solutions which are nearly identical at 
the initial time do not remain nearly identical as time progresses, but 
eventually become as different as two solutions chosen randomly at the same 
time and day of the year. The atmosphere contains some periodic oscilla- 
tions - principally the annual and diurnal variations and their overtones 
- which are predictable at essentially infinite range. However, accurate 
long range prediction of the remaining features is impossible because the 
initial state of the atmosphere will always be imperfectly known. 

Lorenz (1969a) suggested three approaches that might be used to esti- 
mate the predictability of the atmosphere: 

a. The dynamical approach wherein a system of equations closely resem- 
bling the atmospheric equations is integrated from slightly different 
initial conditions and then the rate of amplification of the differences 

is determined; 

b. The empirical approach which utilizes similar weather types, usually 
referred to as "analogues" and determines their subsequent separation in 
statistical terms as a function of time; 

c. The empirical-dynamical approach which utilizes a new set of equa- 
tions (derived from the atmospheric equations) which describes a spectral 
distribution of forecast errors. 



39 



Most of the studies which follow (a) used general circulation models 
which had been integrated until they were in approximate balance. At this 
time another integration was begun with an initial state which consisted 
of the balanced field plus a small departure. As this solution departed 
from the control experiment the growth rate of the error field and the 
range of predictability were obtained. Experiments by Mintz and Arakawa, 
Smagorinsky and Leith have been summarized by Charney, et al (1966). 

Charney estimated the RMS doubling time of small errors to be about 4 to 
5 days and a predictability limit imposed by typical observational errors 
of about two weeks. Smagorinsky (1969) presented the results of experi- 
ments carried out at the NOAA Geophysical Fluid Dynamics Laboratory with 
their 10- level primitive equation general circulation model. A random 
temperature disturbance with a standard deviation of 0.5°C at all levels 
was added to the control run. The vertical average of the standard devia- 
tion dropped from 0.5°C to 0.2°C after one day, reflecting a "geostrophic 
adjustment" between the initial disturbed wind field and the undisturbed 
wind field. Thereafter, the error growth was exponential for the next 
seven days with a doubling time of about 2.5 days. Smagorinsky concluded 
that the deterministic limit of predictability for synoptic scale disturb- 
ances is about three weeks; however, the current practical limit is about 
one week. He also noted that short-wave predictability decays most rapidly. 

The cause of the exponential growth of small errors in these studies 
was attributed to baroclinic instability [Charney, et al (1966)]. However, 
Lorenz (1972) suggested that the error growth might be due to barotropic 
instability of wave disturbances. In particular he investigated the 
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stability of a finite amplitude unbounded Rossby wave. He found instability 
when the waves are sufficiently strong or the wave number is sufficiently 
high, and the growth rates are comparable to the growth rates obtained from 
predictability studies. Lilly (1973) emphasized the importance of baro- 
tropic instability in a predictability study which employed a 2-dimensional 
model. Recently, Lorenz (1973b) has concluded that baroclinic instability 
is the most important cause of lack of predictability in the atmosphere. 

Lorenz (1969c) utilized approach b to obtain the rate of separation of 

two fields which were initially similar. Five years of twice-daily height 

values of the 200- , 500- , and 850-mb surfaces on a grid of 1003 points were 

obtained. A weighted root -mean-square height difference is used as a measure 

h tf 

of the difference between two states, or the error. For each pair of states 
occurring within one month of the same time of year, but in different years, 
the error was computed. Numerous mediocre analogues were found but there 
were no really good ones. The smallest errors had a doubling time of about 
eight days. Since larger errors grow less rapidly this is probably an over- 
estimate of the doubling time. Extrapolation with the aid of a quadratic 
hypothesis indicates that very small errors double in about 2.5 days. This 
compares very well with the numerical results reported by Smagorinsky (1969). 

In a more recent paper Lorenz (1973a) has used the same data set to in- 
vestigate the range of predictability. States of the atmosphere separated 
by 12 days or less are found on the average to resemble each other more 
closely than randomly selected states, even after an adjustment for the 
seasonal trend has been made. Higher correlations were obtained with a 
form of damped persistence. These results demonstrate the existence of 
partial predictability of instantaneous weather patterns at least 12 days 
in advance. 
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Lorenz (1969b) followed approach c with a statistical treatment of the 
2-dimensional vorticity equation. He derived an equation for the "error 
energy" which is obtained from the difference between two solutions. 
Thompson (1957) considered similar equations. The linearized equations 
are written in spectral form and ensemble averages are taken. A further 
assumption is included to close the set of predictive equations. The equa 
tions require the specification of a mean energy spectrum, and Lorenz used 
the -5/3 law for the smaller scales. The equations were integrated numeri 
cally from an initial error distribution; the error energy of each wave 
number was not allowed to exceed the mean atmospheric energy for that 
wave number. The numerical solutions show a rapid growth of the error for 
the very small scales and a very slow growth for large scales. Lorenz's 
solutions show that the range of predictability for cumulus scales is 
almost an hour, while the synoptic-scale motions can be predicted a few 
days ahead. Predictability for the largest scale disappears after about 
17 days. In particular if the initial error is confined to the smallest 
scales the error in those scales rapidly reaches the maximum error. Then 
growth occurs in the next larger scales until the error is propagated to 
the largest scale. Lorenz points out that if the error in the smallest 
scale is reduced then the range of predictability of the largest scale 
features will only be increased by a time interval which is less than the 
range of predictability of the small scales where the error was reduced. 
Although these results are dependent on the closure assumption and the 
energy spectrum which is used the study shows clearly the importance of 
the nonlinear propagation of errors between different scales of motion. 

Robinson (1967, 1971) assumed that the dynamic equations do not allow 
one to predict motions of a given scale over periods longer than the fluid 
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elements of this scale maintain their identity against turbulent diffusion 
by smaller scale motions. Then based upon the dissipation characteristics 
of the atmosphere, he derived predictability ranges of a few days for 
synoptic-scale motions and about one hour for cumulus scale motions, which 
are roughly the same as obtained by Lorenz (1969b). 

Leith (1971) and Leith and Kraichman (1972) have also used homogeneous 
isotropic turbulence models to study atmospheric predictability. Leith 
(1971) considered two dimensions and closed his equations with the eddy- 
damped Markovian approximation. Leith and Kraichman (1972) considered both 
two and three dimensions and they used Kraichman' s (1971) closure method 
which is based on the test-field model. The solutions in three dimensions 
showed that even where there is a strong energy cascade there is also an 
upscale propagation of error. The two-dimensional solutions used a "-3" 
power energy spectrum for the smaller scales. Charney (1971, 1973) showed 
that quasi-geos trophic three-dimensional flow should have a M -3 M power 
energy spectrum under proper conditions. Leith and Kraichman (1972) found 
that the two-dimensional flow was more predictable than three-dimensional 
flow. They concluded that an initial state determined with a horizontal 
resolution feasible with a satellite-based observing system would result 
in significant predictability of large scale motion for more than one week. 
They also point out that the test-field model probably underestimates rather 
than overestimates predictability times. 

In summary, the studies of predictability suggest that specific weather 
patterns and events are predictable, or partially so, only for a period of 
at most several weeks; however, the possibility of predicting general trends 
of perhaps temperature and precipitation above or below normal for longer 
periods is not precluded as yet. 
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12. Stochastic Dynamic Prediction 

Recent studies in stochastic dynamic prediction are closely related to 
atmospheric predictability. Epstein (1969) developed a method for includ- 
ing the effect of random errors in initial conditions in a forecast model. 
The forecast equations which he treated were of the spectral type with 
quadratic nonlinearities. He first took an ensemble average of the basic 
equations which leads to predictive equations for the ensemble average or 
the expected value of each dependent variable. These equations also in- 
clude the variances and covariances between the dependent variables. Equa- 
tions for the time change of the latter quantities which are second moments 
are obtained by multiplying each equation by each of the variables and 
carrying out the ensemble averages. These equations, however, contain 
the third moments. Equations for the time change of the third moments are 
obtained in the same way and they contain the fourth moments. This is 
merely the closure problem of classical turbulence theory. 

Epstein closed his equations by neglecting the third moments about the 
instantaneous mean. He integrated the resulting equations for the Lorenz 
(1960) three component system. Initially he specified the expected value 
of each variable as well as its variance due to possible observational 
errors. The initial covariances were set equal to zero. A deterministic 
forecast was also performed with original prediction equations. Numerical 
integrations showed that the deterministic solution eventually diverges 
from the stochastic prediction of the expected value. The latter should be 
a better forecast in the statistical sense. The predicted variances grow 
in time which gives the influence of the initial error. Also a large number 
of deterministic forecasts were made with slightly different initial condi- 
tions. These numerical solutions were then averaged using a Monte Carlo 
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method. The difference between the Monte Carlo average and the approximate 
stochastic solution is a measure of the error which is caused by the closure 
assumption. In some of the solutions these quantities remained close while 
in others they began to depart after a few days. 

Fleming (1971a) has continued this work with a two-level quasi-geos trophic 
spectral model similar to the one used by Lorenz (1965). For interpretation 
he divided the energy into a "certain" energy and "uncertain" energy. He 
also reconsidered Epstein's closure assumption which neglected the third 
order moments. Fleming tested the quasi-normal approximation which computes 
the fourth order moments in the third order equations in terms of the second 
order moments. This requires the numerical solution of the third order equa- 
tions in addition to the others. He found that this approximation was better 
up to a certain time, after which it became unstable. He then considered a 
modified form in which a linear damping term was added to the third order 
equations. This eddy-damped quasi-normal approximation was stable and gave 
excellent results in most cases. 

Fleming (1971b) used the eddy-damped quasi-normal approximation in a 
stochastic formulation of the barotropic model. He used two wave numbers 
in latitude and 15 in longitude in a study of predictability. The predicta- 
bility times obtained are similar to those obtained by Lorenz (1969b). Pre- 
dictability studies using the baroclinic model of the earlier paper showed 
that predictability is increased when heating and friction are present. 

Fleming (1972) considered the effect of random variations in the thermal 
forcing and the effect of errors resulting from the improper treatment of 
the smaller scales. 



45 



Leith (1973) has concluded that the stochastic dynamic method cannot 
be used for the larger numerical models because of the very excessive com- 
puter time requirements. He suggests instead the use of the Monte Carlo 
procedure which involves a collection of deterministic forecasts. The 
procedure was evaluated with a simple two-dimensional turbulence model of 
large-scale atmospheric motions. A study of the dependence on sample size 
of mean square forecast error, both with and without a final linear regres- 
sion correct ion, showed a considerable improvement with moderate sample 
sizes over conventional single forecasts without regression. These results 
suggest that a sample size of about 10 may be adequate for producing the 
error variance information needed for optimal data assimilation. 
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13. Spectral Methods 



The spectral method expresses the spatial variation of the prediction 
fields in terms of a series of orthogonal functions. The coefficients in 
the series are now the forecast quantities rather than gridpoint values of 
the original dependent variables. We will illustrate the technique with 
the barotropic vorticity equation: 



where t|r is the stream function and the beta term has been neglected for 
simplicity. We expand \|/ into the following series: 



The integration in (21) is carried out over the entire forecast region. 



where the eigenvalues K a are positive and increase for the smaller scale 
eigenfunctions. Substitute (20) into the forecast equation (19), multiply 
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The functions Y a are orthogonal and normalized so that 
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the interaction coefficients, which can be computed once and for all, are 
given by 

L 7,a,|3-- K eS Y 7* lk •™ a * v v s • (24> 

Equation (23) represents an infinite number of ordinary differential equa- 
tions when all appropriate values of 7 are considered. In practice the 
series given by ( 20 ) is truncated in such a way that the desired meteorologi- 
cal features are properly described. The equations for the remaining 
coefficients can then be integrated in time numerically. 

Lorenz (1960) treated equation (19) with this procedure in cartesian 
coordinates. In that case the eigenfunctions are products of sines and 
cosines and the interaction coefficients take on a particularly simple form. 
He noted that when the series is truncated total energy will still be con- . 
served if the excluded coefficients are set to zero in the interaction suras. 
He also demonstrated the usefulness of very low order systems with his study 
of a 3-component set. 

Silberman (1954) first applied this technique in spherical coordinates 
although he kept the zonal mean flow fixed in time. Platzman (1960) and 
Baer and Platzman (1961) performed a complete treatment of the barotropic 
vorticity equation in spherical coordinates and also carried out some 
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experiments. In spherical coordinates the eigenfunctions are spherical 
harmonics which are composed of sines and cosines of the longitude and 
Legendre polynomials of the sine of the latitude. With these functions 
the interaction coefficients are more complicated and they have more non- 
zero elements than in cartesian coordinates. 

The spectral method has several advantages over the gridpoint method. 
First the spectral method computes spatial derivatives exactly so that the 
phase error which occurs with the finite difference method is eliminated. 
Also the aliasing that occurs with finite differences is excluded and as 
a result it is easy to conserve certain quantities which are conserved in 
the continuous equations. Poisson equations are easily solved without 
relaxation or other inversion techniques because of relation (22). Another 
advantage is the treatment of global motions without the presence of 
singularities . 

The most important disadvantage of the spectral method is that it 
requires much more computer time than the gridpoint method, if there are 
very many degrees of freedom. This can be seen from equation (23) which 
shows that for each degree of freedom many products must be computed in 
the nonlinear term. With the gridpoint method, for each degree of freedom, 
there are only a few products involving quantities at surrounding grid- 
points. The spectral method is generally more complicated to apply and 
it is not convenient for complicated boundaries. The spectral method also 
suffers when values must be evaluated locally such as in determining con- 
densation criteria. We shall see that some of these difficulties have 
been alleviated with recent developments. 
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Robert (1966, 1968) modified the spectral technique in two studies with 
the primitive equations on the sphere. In place of the spherical harmonic 
functions he used the following function of sines and cosines: 

G m (A,cp) = e 11 ^ cos ^ cp sin n cp . (25) 

n 

These functions are not orthogonal, but the spherical harmonic functions 
can be written as a sum as of the functions given by (25). The nonlinear 
terms are easily computed with these functions, and the number of elements 
involved is much less than with the spherical harmonics. When the 
orthogonality relations are required the spherical harmonics can be formed 
in terms of these functions. Using this simplified procedure Robert 
carried a number of spherical integrations with a relatively small number 
of terms. 

A very important advance for the spectral method was the development of 

the fast Fourier transform by Cooley and Tukey (1965). Their technique 

allows Fourier transformation of a field with N degrees of freedom with 

2 

order N N operations, while the direct method requires order N 

operations. This allows rapid calculation of local fields by summing 
the series with the fast Fourier transformation, such as might be required 
for condensation tests. 

Orszag (1969, 1970) has used the fast Fourier technique to save compu- 
tation time in the evaluation of the nonlinear terms. He carries out all 
differentiation of quantities while they are represented by the series, 
but products of fields are performed with the fields at gridpoints. This 
requires an inverse Fourier transform to obtain data at gridpoints and a 
regular Fourier transformation to return to spectral form. The fast Fourier 
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transformation can be used in both of these operations with a great time 
savings. When this procedure can be used it greatly reduces prediction 
time, for large numbers of variables, because the sum in (23) is replaced 
by a much faster operation. 

Merilees (1973a) has developed an algorithm for computing the sum of 
a series of spherical harmonics. The algorithm is approximately 10-20 
times faster than the standard method, but it suffers a precision problem 
and it breaks down when the resolution is too high. When this method is 
used^ following Orszag, for rapid computation of the nonlinear terms, the 
prediction time for a global spectral model can be greatly reduced. 

Bourke (1972) formulated a global spectral model based upon the one- 
layer shallow water equations. He used the vorticity and divergence equa- 
tions instead of the 2 components of the equation of motion. The technique 
for evaluating the nonlinear terms which was developed by Orszag (1970) 
and Eliasen, et al (1970) is employed. The computational efficiency of 
the model was found to be far superior to that of an equivalent model based 
on the interaction coefficients. This model, in integrations of 116 days, 
satisfied the principles of conservation of energy, angular momentum and 
square potential vorticity to a high degree. 

Daley (1973) used Bourke's (1972) model to examine the possibility of 
using different spatial resolution for different dependent variables. He 
suggested that for the smaller scales it would be sufficient to predict 
only the wind field because for those scales the pressure field adjusts to 
the wind field. In a test the smaller scales were predicted with a filtered 
model and the larger scales with the full primitive equations. He found 
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that this combined system gave better forecasts than a primitive equation 
system with intermediate resolution. The application of this procedure to 
baroclinic equations may be more difficult. 

Recent developments have made the spectral models much more competi- 
tive with gridpoint models with respect to computational efficiency. For 
the same resolution the spectral models can be expected to give better 
forecasts, and the forecasts may be better even with somewhat poorer reso- 
lution. A series of tests are needed to determine whether or not better 
forecasts can be made with the spectral method with the same computational 
effort. 

Orszag (1972) and Merilees (1973b) have proposed a forecast technique 
which is a combination of the finite difference method and spectral method. 
This technique which is known as the pseudospectral approximation employs 
finite differences except in the computation of spatial derivatives. 

Before the derivatives are computed, the gridpoint data is fast Fourier 
analyzed. Continuous derivatives are then computed and the series is summed 
with the fast Fourier transform. This procedure is used for both longi- 
tudinal and latitudinal derivatives. The latitudinal derivatives in some 
cases require a sign change at the poles. Merilees (1973b) found excellent 
results in a test prediction of Haurwitz waves with the shallow water 
equations . 

14. Finite Element Method 

The extensive success of the finite element method in solid mechanics 
has attracted the interest of investigators in other fields. Salinas (1973) 
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has given a brief review of the method which proceeds as follows: At the 

start, an approximate solution in the form of a linear combination of speci- 
fied (basis) functions is assumed. The coefficients (multipliers) of the 
basis functions are to be determined to yield the best approximate solu- 
tions. This is accomplished by minimizing a measure of the error (called 
the residual function) associated with the assumed solution. Several tech- 
niques are available for minimizing the residuals. Normally the basis 
functions are defined in such a way that they have a simple variation 
(quadratic, cubic, etc.) over an element of area (piecewise polynomials), 
outside of which they are zero. The "elements" can have a variety of 
shapes. The finite element method is a generalization of the method of 
weighted residuals. 

Newton (1973) has successfully applied the technique to gravity waves 
radiating out from an oscillating ship. Gallagher and Chan (1973) have 
treated the steady state circulation in Lake Ontario. Triangular and 
otherwise shaped elements were used to accurately fit the lake’s coastline. 
Wang, et al (1972) applied the method to the one-dimensional shallow water 
equations. They obtained better results for both the linear and nonlinear 
cases than with the usual finite difference methods. Price and MacPherson 
(1973) have applied the technique to the two-level general circulation model 
developed by Mintz and Arakawa (Langlois and Kwok, 1969). They arranged 
the elements in such a way that the area of each element is nearly con- 
stant over the entire globe. They also provided for a subregion in which 
the elements telescoped to a smaller size. Predictions with this model 
gave encouraging results. 

For meteorology, the principal advantage of the method as shown by 
Price and MacPherson (1973), is the possibility of easily changing the 
element size and shape. This would be useful in situations where meshed 
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models are now used and also on the sphere where the elements could be kept 
at a fixed scale. 

15. Summary 

Research on initialization has intensified in the last few years. Most 
operational forecast models employ the primitive equations which are very 
sensitive to the initial balance between the wind and the pressure fields. 
Short range prediction of precipitation is also very sensitive to the initial 
divergence field. The development of global models requires proper initiali- 
zation in the tropics which may differ from middle latitudes. The most 
promising initialization techniques include the variational formulations 
and the averaging of backward and forward forecasts about the initial state. 
The assimilation of nonsynoptic data into a forecast model has received con- 
siderable attention. This has been motivated by availability of nearly con- 
tinuous satellite data. The main problem here involves controlling the 
imbalance introduced into the forecast by the localized new data. 

Several groups have developed limited area models which are designed to 
give better forecasts in specific regions. These models use boundary condi- 
tions from hemispheric or global models and in some cases the limited area 
forecast is allowed to affect the exterior region. Isentropic models have 
been developed and applied to smaller scale synoptic features. 

Operational predictions of pressure and wind have shown continuing im- 
provement as a consequence of progressively better numerical models. How- 
ever, the concomitant improvement in forecasting precipitation and tempera- 
ture has not kept the same pace. It is expected that forecasts of weather 
elements will improve as the limited area models produce better descriptions 
of mesoscale systems. 
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Studies of predictability suggest that specific weather patterns and 
events are predictable, or partially so, only for at most several weeks. 

The possibility of predicting general trends above or below normal for longer 
periods is not as yet precluded. Perhaps the consideration of atmosphere- 
ocean interaction will lead to longer range prediction of weather trends. 

The stochastic dynamical prediction procedure is closely related to 
studies of atmospheric predictability. The procedure uses the initial un- 
certainty in the initial state to predict the uncertainty in the forecast 
at later times as well as the expected value. This method cannot be used 
for operational forecasts in the near future because of the large computer 
time requirements. However, an indication of the growth of uncertainty 
can be obtained by examining a relatively small number of deterministic 
forecasts with slightly different initial conditions. 

The development of the semi-implicit method is a major advance in finite 
differencing. The method treats the gravity wave terms implicitly and the 
advection terms explicitly. This permits the use of a longer time step 
with the additional requirement that an elliptic equation be solved at each 
time step. The net effect of this is a reduction in computer time by a 
factor of at least 2. This has lead to the development of new methods for 
solving elliptic equations which are faster and more accurate than the 
traditional relaxation methods. 

Present operational forecasts often underpredict the movement of synoptic 
disturbances and this error has been attributed to space truncation error 
in the finite difference equations. Several groups are experimenting with 
fourth order space differencing in order to reduce this error. Also the 
pseudospectral method has been put forward as a method to eliminate space 
truncation effects . 



55 



New interest has developed in the use of the spectral method since the 
introduction of the fast Fourier transform. The use of this transform 
greatly speeds up computation time with the spectral method. Further 
tests are required to determine whether the spectral method will give 
better forecasts than the finite difference method for the same amount 
of computer time. 

The finite element method of solid mechanics has been recently intro- 
duced into meteorology. This method has the advantage of a very flexible 
element size, but it has not been widely tested on meteorological problems. 
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A. PRE NWP 

B APRIL 1958 TO FEBRUARY 1964 BASED ON NWP 500-MB BAROTROPIC 
ANO 3-LEVEL BAROCLINIC GUIOANCE 

C. MARCH 1964 TO MAY 1966. BASEO ON REEO SURFACE GUIDANCE. 

D. JUNE 1966 TO AUGUST 1967 BASEO ON P E SURFACE GUIOANCE. 

E SEPTEMBER 1967 TO DECEMBER 1971. BASEO ON P.E. SURFACE GUIOANCE 

□ INOICATES HUMAN IMPROVEMENT. 

Figure 2.--NMC average 30-hour surface scores 




A PRE NWP 

B APRIL 1958 TO JUNE 1962. BASED ON BAROTROPIC GUIDANCE. 

C. JULY 1962 TO JUNE 1964. BASED ON 3-LEVEL BAROCLINIC GUIDANCE. 

D. JULY 1964 TO JANUARY 1971. BASED ON 3-LEVEL BAROCLINIC GUIDANCE 
TO JUNE 1966 AND P.E BAROCLINIC GUIDANCE THEREAFTER. 

E JANUARY 1971 TO DECEMBER 1971. BASED ON P.E. BAROCLINIC GUIDANCE. 

□ INDICATES HUMAN IMPROVEMENT 

Figure 3.--NMC average 36-hour 500-mb scores 
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